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UNIVERSITY OF ILLINOIS 
DIGITAL COMPUTER LABORATORY 
ILLIAC PROGRAM LIBRARY 

Library Routine M 25 - 26g 

Eigenvalues and Eigenvectors of a Symmetric Matrix (SADOI Only) 

Open * 

155 + (Hi) = 164 

0-21 

3 
5n milliseconds per iteration; the number of iterations 

varies from k to 8. 

The sums of the squares of the elements of the matrix must 

be less than one half. 

About 8 or 9 decimal places 

S3: Location of symmetric matrix 

S^: OOF OOnF n = order of matrix. 

S5: OOF OOmF m = 2n if eigenvectors are wanted, otherwise 

m = n. 

S6: OOF n >g +1 ^ S3 

S7: OOF 00n 2 S6 

SYMBOLIC ADDRESSES USED: 

M 25, R 1 

The lower off diagonal elements and diagonal elements should 

be stored consecutively beginning at location S3. If only 

eigenvalues are desired then only n(n+l)/2 memory locations 

are required for the matrix A. If eigenvectors as veil 

are wanted then an additional n 2 memory locations must be 

reserved for the eigenvectors. 

The off diagonal elements of the original matrix are 

reduced to close to zero. (See Description of Method 

and Convergence Criterion). The elements of the 

orthogonal matrix of eigenvectors are stored consecutively 

in locations beginning at S6, scaled by l/2. Thus the 

matrix 

a ll 
8 21 a 22 



METHOD OF USE: 
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nl 



nn 



becomes 
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e l 
e 2 

e.. 
3 



.... e followed by 
II 



k ll' k 12' '"'*]* 
k gl , k 22 , ..., k^ 

k .......... *_ where # we denote by v. th« Jth 

eigenvector its n components are: k,,, ^p^ •••# * <> 
i.e. each column is an eigenvector. 
MOTES: (l) If there is an arithmetic error (see below) two 

■F's will be punched out and the machine will stop on 
an OF at the right hand side of the 158th word of this 
routine. The computation can be continued by raising 
and lowering the white switch. 

(2) If the original matrix was scaled down by k, then 
the eigenvalues are scaled by k. 

(3) R 1 is at location 156L. 

(k) This program is a revision of M - l4l and replaces 
it in the library. 

(5) Care should be taken so that the user does not exceed 
the capacity of the Williams Memory. In the complete 
program version of this routine the maximum size of the 
matrix A is 21 x 21 if eigenvalues and eigenvectors are 
desired and 37 * 37 if eigenvalues only are desired. If 
eigenvalues and eigenvectors are desired then the Williams 
Memory from S3 to S3 + n(n+l)/2 + n - 1 will be reserved 
for the matrix and the eigenvectors. If eigenvalues only 
are desired then Williams Memory S3, S3 + n(n+l)/2 - 1 
will be reserved far the matrix A. 

(6) Since IEa.jjb.jJ < E|a ± . | E|b,jJ *= 1 the results of 
the successive multiplication of orthogonal matrices 
remain in scale if round-off is ignored . In the presence 
of round-off, overflow is prevented by using l/2 I (the 
identity matrix), instead of I. 
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BRIEF DESCRIPTION OF THE METHOD: 

Let A be a symmetric n x n matrix. We define a matrix 

as follows: 

k j 

v where all diagonal elements 

except 0.., VV equal 1 and 

°<5 k * A ... ^^''•r.nflf / a11 off diagonal elements 

except and equal 0. 

It is easily verified that , is an orthogonal matrix and 

m i J* ■>, 

hence 0., = 0., . Since A is symmetric then, by a well • 
jk jk 

known theorem of algebra, there exists an orthogonal matrix 

T 
such that AO = D a diagonal matrix; in the language 

of linear algebra "any real quadratic form in n variables 
assumes the diagonal form, relative to a suitable ortho- 
normal basis." (Birkhoff and Maclane p. 277)- This routine 

T 
constructs the matrix by an iterative process which at 

the same time reduces A to its diagonal form D. In this 

iterative procedure the orthogonal matrices , play a 

fundamental role. 

T 
First we form 0, v AO., = B and note that B differs from 
J K Jk 

A only in the jth and kth rows and jth and kth columns. 

a. , cos "f + a. sin | 
fb.- = - a. t sin I + a. cos i 

^kk = a jJ co£ * ^ + 2a jk cos ' sin ' + a kk sin T 

[b = a + a,, - b,, [under an orthogonal transformation 
* jj jj kx xk 

frjk = a jk^ cos2v ^ " s±n2 ^y + ( a jj " a kk^ sin t c ° s f 

[ = a cos 2^ + 1/2 (a ^ - a^sin 2f 

A necessary and sufficient condition for b,, = is that 

"T satisfies the following equation: 

2a., 

ten2 f = ^J* 

a tt " a JJ 

It should be noted that T is not uniquely determined by 
this condition as T + ™/2 will also satisfy this equation. 
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But the important thing is to choose sin \ and cos *f 
in a manner consistent with the condition 

tan2t- ^ 

and no matter how this is done the result is to set 
b.,= 0. We have shown that given any non-zero off 
diagonal element a., an orthogonal transformation 
0., may be constructed auch that in the matrix CLjJIO - B, 
then b„ =0. If we choose 0.. in this way it is 

easily verified that 

r» 2 v 2 2 

^ b rs = £• a - 2a i.e. the sums of the squares 

r^s r/s J 

of the off diagonal elements are reduced by a positive 

2 
amount 2a,,. This is easily verified as follows: 
J^ 2 2 2 ? 

for /fj, k bw,+b^ = a kJ? + a j/ " Therefore > 

£ b 2 = E a 2 - 2a 2 
r£ rS r/s rs ^' 

The method is therefore to select in some way a sequence 

of a , / (j f k) and reduce these in succession to 

zero via the method described above. After a finite 

number of these rotations the off diagonal elements will 

be sufficiently small (zero inside the machine since the 

sura of their squares is constantly being reduced by a 

positive amount) and the process has converged. Let D 

represent the matrix in the machine to which A has been 

reduced and let D be the matrix obtained from D by setting 

all off diagonal elements to zero. Then if we denote by 

X, the eigenvalues of D and by X, the eigenvalues of D 

then by a well known theorem of Courant |x, - X, | < |D - D| 

i — i J J 

where |D - D| refers to the matrix norm. 

If v is an eigenvector with respect to D, its representation 

with respect to A must "be found. A0v = Xv implies 

(multiplying on the left by 0) that A(0v) = X(0v)j hence 

if v is an eigenvector with respect to D, Ov is an 

eigenvector with respect to A. Since (l, 0, ..., 0), 

(0, 1, 0, ..., 0), ..., (0, ..., 0, l) are the eigenvectors 
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with respect to D, It follows that the columns of the 
Matrix represent the eigenvectors with respect to A. 
Thus if we want the eigenvectors all we need do is 
form T i+1 = ^B. (where T = I) after each iteration and 

B T ... B-f A B, B ... B = D. 

m i. i. 2 m 

We note in passing that the passage from D to D is not 
without its dangers as the following simple example shows. 

- * ■(>;?) —GO 

then the sum of the squares of the off diagonal elements 
of A is e/2 hut the eigenvectors of A are not 

detailed discussion of these and other matters relating 
to this method see the paper by H.H. Goldstine, F.J. 
Murray, and J. Von Neumann: The Jacob! Method for Real 
Symmetric Matrices , Journal of the ACM, volume 6, 
January, 1959 j no. 1. 
DESCRIPTION OF NUMERICAL OPERATIONS PERFORMED; 

When the routine begins, the original symmetric matrix 
A = A^ is stored in the consecutive memory locations, 
S3, S3 + 1, •••, etc. and if the eigenvectors are desired 
the unit matrix, T Q , is generated and stored in the 
Williams Memory beginning at S6 and ending at S6 + n - 1. 
After executing the ith transformation each element of 
A . , occupies the memory locations previously occupied 
by the corresponding element of A. and the same applies 
to T. , and T. . Actually, the machine changes only those 
elements affected by the matrix multiplication and leaves 
the remaining elements unaltered. The ith transformation, 
reducing element a., to zero, will alter only the elements 
of A. in the jth and kth columns and rows and will alter 
only the elements of T. in the jth and kth columns. This 
means we start with a pair of elements, a, .. , a^, and 
work our way across the kth and jth rows until we reach 
the diagonal and then go down the kth and jth columns until we 
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reach the last row of the matrix T (if eigenvectors are 
desired; otherwise, we go until the last row of the matrix 
A. is reached). In order to keep track of our progress 
we use jf as a tally and when /f - S5 is positive we are 
through. The elements are transformed according to the 
equations 

w *fc +1) -2}°** + *$'^ 

where a* ,, = a^ j a * = a ^ . However, there are three 

special cases which do not use these equations and since 

these elements have already been computed incorrectly we 

merely write the correct values over the incorrect* ones. 

We set a = and transform the diagonal elements a., 
jk JJ 

and a,, according to the equations 

<5> a £ +1) - ^ ^ + 4 } si » 2f + a S } sln2 * 

I 4 ' a jj " a kk + a jj a kk 

The elements a., for successive transformations are selected 

Jk 
in consecutive order along successive rows, i.e., we use 

a in matrix iL , a.., in matrix A , etc. Of course, if 

a., = already we do not rotate but go on to a j* k:+1 ' We 

define one iteration to be the (n - n)/2 transformations 

required to reduce each off diagonal element to zero once. 

It should be stated that once an element is reduced to 

zero it will not> in general, remain zero during subsequent 

transformations. However, the sum of the squares of the 

off diagonal elements will be decreased each time by an 

amount equal to 2a. v and thus will become small (zero 

inside the machine), after a finite number of transformations. 

The program can be divided roughly into two parts: 

(l) Computation of sin"f and cos^ I 

We always choose cos V P> 0. Define a = 2a. k ; b = a^ - a.y 

Case 1: |a | < |b | then m = tan 2^f 

m _ a 
2 " 2b 
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q - 1 








| cos 2f = 


S 


sin l cos [ 


« ^ cos^ 




Case 2: |a | > 


|b| then m = 


cot 2^ 




m b 

2 " 2a 


q m 






1 sO 


Isl 






^ cos 2 J = 


rr + ^ 






Case 2a: S < 0, 


1 
r = "IT 






Case 2b: S > 0, 


1 

r = z 








sin"Y cos \ 


r 






[1 m2 





In both cases cos T = / s + p cos 2T 

sin"^ = [sin V cos f ]/cos / 

(2) Matrix multiplication: 

This has already been mentioned. However 3 we need to 

describe the method of obtaining the addresses of the 

various elements involved in the matrix multiplication. 

In order to simplify the discussion we make a slight change 

of notation - we begin our count from instead of 1; i.e., 

we denote an element of the first row by a * instead of 

a, g and similarly an element from the first column by 

a n „ instead of a n . 

First of all once a., is chosen we need to compute (kj 0)j 

J* 
(k, k); (j, 0); (j, j) where these denote the addresses 

of a 1A ; a. : a.„: a. , . 
k0* kk* JO' jj 

(k, 0) - | (k 2 + k) + S3 

(k, k) =■- (k, 0) + k 

(J, 0) = | (j 2 + J) + S3 

Us J) = Us 0) + j. 

These addresses are planted in the orders which carry out 

the operations indicated in equations (l), (2), (3)j (*0 

above. How (k,/) and (j, J( ) j(= 0, 1, ..., S5 should 
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be increased by 1 as one moves across the kth and jth rows 
mrtil tn e diagonal elements are reached. Then they should 
be increased by / as one moves down the jth and kth columns 
until the last row of A is reached, i.e., J? . a - l. 
Then if eigenvalues are desired they should be increased 
by n until tf = S5 . We handle this by storing two 
increments in the address portions of 150 (M 25). The 
increment in the left hand address is used to move along 
the path starting at (k, 0) and the increment in the right 
hand address is used to move along the path starting at 
(j, 0). The increments are not always the same since one 
path reaches the diagonal sooner than the other. The 
determination of the increments requires that our tally 
A (in location 11F) be compared with k and j so as to change 
the increment from 1 to^ along each path. It is also 
necessary to compare / with n so as to change the increments 
from J to n. Finally we compare ^with S5 in order to 
get out of this loop. As was mentioned earlier we take 

care of the three special cases (a^ i+1 ' . a^ i+1 ^ fl( i+1 ) \ 

v Jj ' kk > a jk ' 
after we leave the loop. 



DESCRIPTION OF CONVERGENCE TEST AND ARITHMETIC TEST: 

If we define 
n 



H 2 = Z a 2 E . £ a 2 S = Z a 2 

then we have 

S - N 2 - (|) E. (1) 

Under an orthogonal transformation N 2 remains invariant and 
E is reduced by an amount equal to twice the square of the 
element which goes to zero under the transformation, i.e. 

E i + i = E i- 2ai W 

Thus using (l), 



(N 2 - \ \) + a 



2 
2 



S i +a j'k' (2) 



- 9 - 

and we see that ] S ? forms a monotone-Increasing sequence 
which approaches N 2 as E approaches zero. This is the basis 
of our convergence test. 

Since we form S at the end of each iteration (defined as 
i (n - n) orthogonal transformations) then 

s< i+l) = s^ + w (i+l) (5) 

where IT* ' means the sum of the squares of the elements 
which are reduced to zero "by each of the (1/2) (n - n) 
transformations. Equation (3) gives us our convergence 
test and also provides the means to test the accuracy of 
our computation. We test the quantity S* - S .If 
it is negative we know that the process has not converged 
and we then look at the quantity 

ft- IB* 1 * - S< 1+1) + w( i+1 >| 

where M is our tolerance. If this result is negative the 
machine stops on an OF order from the right hand side of 
138L. It should he mentioned that these quantities have 
been computed using double precision. When S - S 
becomes positive the process has converged. 
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lgr 



LOCATION 


ORDER 


NOTES PAGE 1 




00 K(M25) 









Hi 17F 
Hi 18F 






1 


Hi 19F 








L5 15(M25) 




Compute eigenvectors? 


2 


LO 10L 




> Noi 




36 23L 




< Yesi 




22 ;L 








Hi s6 




Begin clearing 


H 


F5 3L 





it memory locations 




HO 3L 




for 1/2 unit matrix 


5 


LO 11L 








32 3L 




> Not done 


6 


L5 20L 








Ho S6 




Put 1/2 down the 


7 


F5 6L 




diagonal, thus 




LH 12L 




get 1/2 I 


8 


HO 6L 
LO 13L 






9 


36 23L 
26 6L 






10 


00 S5 
00 S5 






11 


K2 3L 
Hi S7 






12 


00 F 

00 sH 






1? 


L5 ^OL 

40 ;.V7 






14 


00 IF 
00 IF 






15 


00 SH 
00 SH 






16 


00 S3 
00 S3 






17 


80 S5 
00 S5 







M 25 



LOCATION 


ORDER 


NOTES PAGE 2 


18 


20 F 
00 F 






19 


00 F 
.00 F 






20 


ho F 
00 F 






21 


00 F 
00 F 






22 


JO S6 
74 S6 








00 K 







F5 19(M25) 
40 19(M25) ' 




Advance iterations 
counter 


1 


4l 20F 
4l 2 IF 






2 


4l 4f 




Set j = 


i 


4l 5F 
L5 ^4-F' 




Set k = (from 99) 


k 


IA l4(M25) 

4o 4f 

LO 15(M25) 




"""'■■ 


5 


32 100L 




J 




50 5F 


from 99' 




6 


L5 5F 






j 


Jk 5F 






7 


00 58F 

42 21(M25) 




>(k, 0) =| (k 2 + k) + S3 


8 


00 20F 








46 21(^5) 






9 


L5 21(M25) 








L4 16(M25) 






10 


46 150L 






li 


L4 5F 
42 29L 




(k, 0) + k = (k, k) 


. i 


46 94L 




1 



M 25 



LOCATION 


ORDER 




NOTES page 


3 


12 


k2 95L 
50 4f 




' 




13 


L5 kF 
Jk kF 




/ 




14 


00 38F 

k2 21(M25) 




W 0) =| (j 2 + J) + S3 




15 


00 20F 

46 21(M25) 




( 




16 


L5 21(M25) 
L4 16(JG5) 








17 


k2 130L 




/ 




18 


L4 4f 
46 96L 
42 28L 




(j, j) = (j, 0) + j 




19 


LO 4F 
L4 5F 


form (J, k) 






20 


42 21L 
42 22L 








21 


42 96L 








22 


L3 F 
56 97L 
L5 F 


by 20 
by 20' 


V ■ 0? 




23 


40 7F 
L5 21F 








2k 


50 7F 










Ik 7F 








25 


L4 20F 
40 20F 








26 


S5 F 
40 21F 








27 


L5 7F 
00 IF 








28 


40 7F 
L5 F 


by 18' 


2a j k - a 




29 


40 9F 




a 






L5 F 


by 11 

- 







M 25 



LOCATION 



30 
31 

52 



36 
'7 



40 



k- t 



ORDER 



ko 

LO 

ko 

L7 
L2 

50 
36 

L5 
66 

S5 
10 

40 

T c: 

26 

L5 

66 

35 

10 

ko 

26 
40 
-0 

75 
LA 
•40 



8F 

9F 

10F 

7F 

10F 

19F 

37L 

7F 

10F 

F 

IF 

11F 

18(M?5) 

40L 

10F 

7F 

F 

IF 

11F 

12 4L 

12F 

11F 

11 F 

18(M25) 

IF 

F 

F 



..»2 (Rl) 

40 14F 

L v -. 12F 

50 10F 

66 14F' 

ni f 

22 llbL 



NOTES 
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M 25 



a k* " a jj 



a - b 



= b 



> Use cot 2 I 5 

< Use tan 2 t 




in 



Ik F 



LOCATION 


ORDER 


ROTES PAGE 5 


kQ 


kO IF 

4l F 




J 2 + | cos 2f = cos i 3 


^9 


22 k9L 
50 i*9L 






50 


22 (Rl) 




cos"f in 2F 




L5 12F 






51 


L0 18(M25) 








32 56L 






52 


L5 12F 
36 5^L 






53 


LI 18(M25) 








22 5^L 






5k 


L5 18(M25) 








50 19F 






55 


66 ikT 
S5 F 






56 


22 58L 








19 IF 


■ 




57 


50 IIF 
74 13F 






58 


00 IF 








kO 1^4-F 






59 


66 2F 








35 F 






6o 


kO 3F 




sin f in 3F 




kl IIF 




Clear tally counter 


6i 


27 70L 
L5 IIF 






62 


LO 15(^5) 
36 100L 






63 


iM- 15(^5) 
ko 131L 






61* 


L5 5F 
LO IIF 






65 


32 67L 








L5 ^F 







M 25 



LOCATION 


ORDER 




NOTES 
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66 


LO 11F 
32 68L 








61 


22 69L 

L5 1^(H25) 








68 


46 13IL 
L5 l4(M25) 








69 


42 131L 
L5 131L 








TO 


L4 130L 
40 130L 








71 


k6 74l 
46 78L 




Plant (k,/) 




72 


46 83L 
42 76L 








73 


42 80L 
42 8lL 




Plant (J,/) 




74 


50 F 


by 7i 






75 


7J 2F 
40 lOF 
S5 F 




a *j? C ° sf +8 « 


Bluf 


76 


50 3F 
74 F 


by 72' 






77 


L4 lOF 
40 F 








78 


50 F 
79 3F 


by 71' 


. 




79 


40 lOF 
S5 F 




a j^ cosr " V 


stnf 


80 


50 2F 
74 F 


*y 73 






81 


L4 lOF 
40 F 


■by 73' 






62 


22 82L 










IF) F 








83 


40 F 


by 72 






J 


\Ij 11F 









M 25 



LOCATION 


ORDER 
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8k 


Ih 1^(M25) 
kO 11F 




Advance A 


85 


LO 17(M25) 




A - S5? 




32 6lL 




> No 


86 


LJ 13F 
kO F 




< Compute 


87 


L9 13F 
kO IF 




(i+1) . a (i+D 


88 


50 7F 




and set a;. =0 
jk 




7J l^F 




89 


kO 10F 
S5 F 






90 


50 8F 
74 F 




f 


91 


lA 10F 
kO 10F 






92 


S5 F 
50 9F 






93 


7^ IF 
IX 10F 






9^ 


ko F 
L5 8f 


by 11' 




95 


IA 9F 








LO F 


by 12 




96 


kO F 


by 18 






kl F 


by 21 




97 


L5 5F 

IA 1MM25) 




Advance k 


98 


ko 5F 








LO ^F 




k = J? 


99 


32 2L 
22 5L 






100 


27 63L 








kl 15F 


from 5 




101 


ki i6f 








L5 l6(M25) 







M 25 



LOCATION 



102 
103 
104 
105 
106 
107 
108 
109 
110 
111 
112 

113 
114 

115 
116 

117 
118 

119 



ORDER 



42 104L 
46 104L 
22 103L 
L5 16F 
50 F 
74 F 
L4 15F 
40 15F 
S5 F 
40 16F 
L5 104L 
L4 l4(M25) 
40 104L 
LO 22(M25) 
32 103L 
L5 18F 
LO l6F 
10 39F 
L4 17F 
LO 15F 
36 9(R1) 
L5 128L 
26 126L 
19 34F 
12. F 
36 ll6L 
92 904F 
OF F 
L5 15F 
40 17F 
L5 16F 
40 18F 
26 L 
L5 13F 
LO 123L 
36 121L 



NOTES 
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by 102 



Form Z 



2 . q (i+l) 



M 25 



g(i) . g(i+l) 

If > then 
process has 
converged 

Arithmetic 
fail stop 



^ 



Bypass arithmetic 
test on 1st 
iteration 



LOCATION 


ORDER 


NOTES 
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120 


LJ 13F 
26 48L 




/Overflow test 




121 


LJ 123L 
40 2F 




/ 




122 


22 50L 








00 F 








123 


3L 4095F 




) 






LL 4095F 






124 


F5 11F 
10 IF 


from 39' 






125 


26 40L 
00 F 








126 


40 112L 
L5 129L 


from 113 






127 


40 113L 
26 116L 








128 


36 9(R1) 
L4 20F 








129 


40 F 
19 34F 








130 


00 F 
00 F 








131 


00 F 
00 F 

(Rl) OOK 




Square Root Routine 




_J 












M 25 



